package scu.maqiang.electromagnetics.homogenization;

import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES2T31;
import scu.maqiang.mesh.*;
import scu.maqiang.numeric.*;

public class ScatteringFinePEC {
    public static void main(String[] args) {
        Mesh2T3 mesh = new Mesh2T3().square2D(384, 384, T3Type.Left);
        mesh.scale(1.5, 1.5);
        mesh.translate(-0.25, -0.25);
        //mesh.toTecplot("FineMesh.dat");

        int N = 8;
        double eps = 1.0 / N;
        int DielDomainLabel = 20;
        ScalarFunc dielDomainFunc = (xy, label, param) -> {
            if(Math.abs(xy[0] - 0.5) < 0.5 && Math.abs(xy[1] - 0.5) < 0.5) {
                double cellX = xy[0] * N - (int)(xy[0] * N);
                double cellY = xy[1] * N - (int)(xy[1] * N);
                return Math.abs(cellX - 0.5) < 0.25 && Math.abs(cellY - 0.5) < 0.25? 1.0:-1.0;
            }
            return -1.0;
        };
        int outerBoundaryLabel = 0;
        ScalarFunc outerBoundaryFunc = (xy, label, param) -> 1.0;
        mesh.setBoundaryLabel(outerBoundaryFunc, null, outerBoundaryLabel);

        int innerBoundaryLabel = 10;
        mesh.removeElements(dielDomainFunc, null, innerBoundaryLabel);
        System.out.println(mesh.extractBoundaryNodes(innerBoundaryLabel).size());
//        mesh.displayBoundarys();

        int PMLDomainLabel = 10;
        ScalarFunc PMLDomainFunc = (xy, label, param) -> Math.abs(xy[0] - 0.5) < 0.625 && Math.abs(xy[1] - 0.5) < 0.625?-1.0:1.0;
        mesh.setDomainLabel(PMLDomainFunc, null, PMLDomainLabel);

        int interfaceBoundaryLabel = 6;
        ScalarFunc interBoundaryFunc = (xy, label, param) -> {
            if(Math.abs(Math.abs(xy[0] - 0.5) - 0.625) < 1.0e-6) {
                return Math.abs(xy[1] - 0.5) < 0.625?1.0:-1.0;
            }
            if(Math.abs(Math.abs(xy[1] - 0.5) - 0.625) < 1.0e-8) {
                return Math.abs(xy[0] - 0.5) < 0.625?1.0:-1.0;
            }
            return -1.0;
        };
        mesh.addInterfaceBoundary(interBoundaryFunc, null, interfaceBoundaryLabel);
        System.out.println(mesh.extractBoundaryNodes(interfaceBoundaryLabel).size());
        System.out.println(mesh.extractBoundaryNodes(innerBoundaryLabel).size());
        mesh.toTecplot("PECHoleDomain.dat");
        mesh.saveMesh(MeshFormat.msh, "HoleDomain.msh");


        //参数设置
        double c0  = 3*1e8;           //m/sec, velocity of light in free space
        double nu0 = 120*Math.PI;          // ohm, intrinsic impedance of the free space
        double e0  = (1e-9)/(36*Math.PI);  // F/m, permittivity of free space
        double mu0 = 4*Math.PI*1e-7;       // H/m, permeability of free space

        // Input parameters
        //******************
        double freq = 300;          // MHz, frequency
        freq *= 1e6;     // Hz, frequency
        double lambda0 = c0/freq;      // meter, wavelength
        double k0  = 2*Math.PI/lambda0; // 1/meter, wavenumber
        double omg  = 2*Math.PI*freq;    // rad/sec, radial frequency

        double mur = 1; // relative permittivity of object (used if scat_type='diel')
        double er = 4;  // dielectric constant of dielectric object (used if scat_type='diel')

        //入射角
        double phii = 45;            // degree, angle of incident field
        double varphi = phii * Math.PI/180;  // radian, angle of incident field
        //入射波函数的负值
        ScalarCFunc incidentWaveFunc = (xy, label, param) -> {
            Complex result = new Complex(0.0, 0.0);
            result.real =  -Math.cos(k0 * (xy[0] * Math.cos(varphi) + xy[1] * Math.sin(varphi)));
            result.image = -Math.sin(k0 * (xy[0] * Math.cos(varphi) + xy[1] * Math.sin(varphi)));
            return result;
        };

        FES2T31 fs = new FES2T31(mesh);
        SRMatrix A = new SRMatrix(fs.GetNdof());
        SCMatrix A_PML = new SCMatrix(fs.GetNdof());

        //在自由空间组装整体矩阵
        fs.assembleStiff(new double[]{1.0/mur}, BVPType.COMMON, A, 0);
        fs.assembleMass(new double[]{-k0 * k0 * 1}, BVPType.CONSIST_MASS, A, 0);

        //在PML区域组装PML矩阵
        fs.assemblePMLStiffMass(mur, k0, -k0 * k0 * 1, BVPType.COMMON, BVPType.CONSIST_MASS,
                A_PML, interfaceBoundaryLabel, outerBoundaryLabel, PMLDomainLabel);
        //整体矩阵
        SCMatrix K = A_PML.add(1.0, A);

        //右端项
        Complex[] RHS = Complex.zeros(fs.GetNdof());

        //整个外边界设为0值
        fs.applyBC_MBN(K, Direct.All, outerBoundaryLabel);
        fs.applyBC_MBN(RHS, Direct.All, new Complex(0.0, 0.0), outerBoundaryLabel);

        //PEC边界，结构区域内部边界
        fs.applyBC_MBN(K, Direct.All, innerBoundaryLabel);
        fs.applyBC_MBN(RHS, Direct.All, new ScalarCFunc[]{incidentWaveFunc}, innerBoundaryLabel);


        DirectSSolver solver = new DirectSSolver(K);
        Complex[] x = solver.complexSolve(RHS, 1, 1.0e-6);

        //入射场负值
        Complex[] uint = mesh.functionInMesh(incidentWaveFunc, null);

        //输出结果, 分别为散射场实数, 虚部, 模, 整体场实数, 虚部, 模
        double[][] result = new double[6][];
        result[0] = Complex.Real(x);
        result[1] = Complex.Image(x);
        result[2] = Complex.abs(x);
        Complex[] utot = Complex.add(x, -1.0, uint);
        result[3] = Complex.Real(utot);
        result[4] = Complex.Image(utot);
        result[5] = Complex.abs(utot);
        mesh.toTecplot("Scattering_FinePEC16.dat", result);

    }
}
